clear;close all;clc;
addpath('fun')

setpar

% =========================================================================
% Cases
scenario = 'FinE';   %(NoFin,Fin,FinE,NoFriction)

switch scenario 
    case 'NoFriction'
fmat_scenario  = 0;
fprob_scenario = 1;
fini_scenario  = zeros(nhh,1);
Vini = -22*ones(np,na,nz,np);

    case 'NoFin'
fmat_scenario  = fmat;
fprob_scenario = fprob;
fini_scenario  = fini;
Vini = -22*ones(np,na,nz,np,nf);

    case 'Fin'
fmat_scenario  = qfmat;
fprob_scenario = qfprob;
fini_scenario  = qfini;
Vini = -22*ones(np,na,nz,np,nqf);

    case 'FinE'
fmat_scenario  = qfemat;
fprob_scenario = qfeprob;
fini_scenario  = qfeini;
Vini = -22*ones(np,na,nz,np,nqfe);
end


% =========================================================================
% I. Solve value functions w. fintech: not depend on hh ini condition
% =========================================================================

V_bf = V(globvar,par,na1,rb_bf,Vini,fmat_scenario,fprob_scenario); % before

Vini = V_bf;
V_af = V(globvar,par,na1,rb_af,Vini,fmat_scenario,fprob_scenario); % after


% =========================================================================
% II. IRF
% =========================================================================

irf_adj_yp  = zeros(nhh,T);
irf_refi_yp = zeros(nhh,T);
irf_a1_yp   = zeros(nhh,T);
irf_b1_yp   = zeros(nhh,T);
irf_c_yp    = zeros(nhh,T);

irf_adj_ypr  = zeros(nhh,T);
irf_refi_ypr = zeros(nhh,T);
irf_a1_ypr   = zeros(nhh,T);
irf_b1_ypr   = zeros(nhh,T);
irf_c_ypr    = zeros(nhh,T);

nw=24
parpool(nw)

% 1/2: yp shock only 
aini = AINI;
bini = BINI;
zini = zini_M;
pini = median(pmat);

for t = 1:T
    disp(t)
    
    parfor i = 1:nhh  
    tic
[adj(i,1),refi(i,1),a1(i,1),b1(i,1),c(i,1)] =...
    getimpact(globvar,par,na1_opt,rb_bf,rb_bf,V_bf,V_bf,bini(i),aini(i),zini(i,1),pini,fini_scenario(i),fmat_scenario,fprob_scenario,t); 
    toc
    end
    
irf_adj_yp(:,t)  = adj;
irf_refi_yp(:,t) = refi; 
irf_a1_yp(:,t)   = a1;
irf_b1_yp(:,t)   = b1;
irf_c_yp(:,t)    = c;

aini = a1;
bini = b1;
end

% 2/2: ypr shock
aini = AINI;
bini = BINI;
RB_BF = rb_bf*ones(nhh,1);

for t = 1:T
    
    parfor i = 1:nhh  
    tic
[adj(i,1),refi(i,1),a1(i,1),b1(i,1),c(i,1)] =...
    getimpact(globvar,par,na1_opt,RB_BF(i),rb_af,V_bf,V_af,bini(i),aini(i),zini(i,1),pini,fini_scenario(i),fmat_scenario,fprob_scenario,t);  
    toc
    end
    
irf_adj_ypr(:,t)  = adj;
irf_refi_ypr(:,t) = refi; 
irf_a1_ypr(:,t)   = a1;
irf_b1_ypr(:,t)   = b1;
irf_c_ypr(:,t)    = c;

aini  = a1;
bini  = b1;
RB_BF = adj*rb_af+(1-adj).*RB_BF;
end


% save file
if T==1 & theta==0.2 % (1) baseline

    switch scenario 
    case 'NoFriction'
filename =['NoFriction','_fH0_x',num2str(-x*100),'.mat'];

    case 'NoFin'
filename =['NoFin','_pf',num2str(pf*100),'_x',num2str(-x*100),'.mat'];

    case 'Fin'
filename =['Fin','_pq',num2str(pq*100),'_q',num2str(q*100),'_pf',num2str(pf*100),'_x',num2str(-x*100),'.mat'];

    case 'FinE'
filename =['FinE','_pq',num2str(pq*100),'_q',num2str(q*100),'_pe',num2str(pe*100),'_e',num2str(e*100),'_pf',num2str(pf*100),'_x',num2str(-x*100),'.mat'];
    end


elseif T==1 & theta~=0.2  % (2) for changing theta

    switch scenario 
    case 'NoFin'
filename =['NoFin_theta',num2str(theta*100),'.mat'];

    case 'Fin'
filename =['Fin_theta',num2str(theta*100),'.mat'];

    case 'FinE'
filename =['FinE_theta',num2str(theta*100),'.mat'];
    end


else     % (3) T>1 for IRF

    switch scenario 
    case 'NoFin'
filename =['NoFin','_IRF',num2str(T),'_pf',num2str(pf*100),'_fH',num2str(-fH*100),'_x',num2str(-x*100),'.mat'];

    case 'Fin'
filename =['Fin','_IRF',num2str(T),'_pq',num2str(pq*100),'_q',num2str(q*100),'.mat'];

    case 'FinE'
filename =['FinE','_IRF',num2str(T),'_pq',num2str(pq*100),'_q',num2str(q*100),'_pe',num2str(pe*100),'_e',num2str(e*100),'.mat'];
    end

end

save(filename);

delete(gcp('nocreate'))

